Equation of state of superfiuid neutron matter and the calculation of 1 S'o pairing gap. 
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We present a Quantum Monte Carlo study of the zero temperature equation of state of neutron 
matter and the computation of the 1 So pairing gap in the low-density regime with p < 0.04 fm~ 3 . 
The system is described by a non-relativistic nuclear Hamiltonian including both two- and three- 
nucleon interactions of the Argonne and Urbana type. This model interaction provides very accurate 
results in the calculation of the binding energy of light nuclei. A suppression of the gap with respect 
to the pure BCS theory is found, but sensibly weaker than in other works that attempt to include 
polarization effects in an approximate way. 



Superfluidity of neutron matter is currently a subject 
of great interest across the astrophysics, nuclear physics, 
and many-body physics communities. Pairing of nucle- 
ons may occur in different channels. At densities rele- 
vant for the phenomenology of the stellar outer core, a 
crucial role is played by the energy gap of the 1 So paired 
state. It has direct consequences on the low density equa- 
tion of state, the cooling process, the phenomenology of 
glitches 0, [3, and on neutrino emission. In addition, the 
effect of pairing of neutrons in the low density regime play 
an important role also in the pairing effects observed in 
neutron-rich nuclei, where the energy is sensibly more 
stable if the total number of nucleons A is even[3j (see 
for example Ref. [H for a review). 

In neutron matter superfluidity of the BCS type occurs 
up to densities of order 2 x 10~ 2 fm~ 3 , corresponding to 
a value of the parameter £ = kpa < —14.8, where a is 
the scattering length, which for a pure neutron matter is 
large and negative (a ~ —18.5 fm). Such value of £ gives 
rise to a weak BCS superfiuid and, at the same time, 
induces strong correlations among neutrons. 

This fact leads to a second main reason of interest, i.e. 
the study of the interplay of the correlations induced by 
the strong repulsion among neutrons at short distance, 
and by the spin dependent forces with a strong tensor 
component, and the occurrence of the BCS state. A clean 
signature of this effect can be found in the behavior of 
the BCS energy gap A. This quantity has been computed 
by several authors within many different approximation 
schemes. A mean-field BCS treatment gives a peak value 
of the energy gap A ~ 3 MeV at &f ~ 0.8 fm -1 , al- 
most irrespective of the NN potential^, 0|. This is not 
surprising because all NN potentials are fitted to repro- 
duce the same S- and P-wave components in the scatter- 
ing experiments. The situation becomes more intricate 
when so-called polarization effects, i.e. the interaction 
with the surrounding medium are introduced. There are 
two ways to attack the problem. The first is still based 
on the solution of the BCS equation(a two body prob- 



lem) with an effective interaction which approximates the 
background effects. Alternatively, a more rigorous calcu- 
lation should include many-body effects by solving the 
many-body problem with the full interaction, and com- 
pute the gap between the BCS and normal state directly 
as the energy difference. At this level, an accurate com- 
putation of the BCS gap can be obtained only if 1) A 
realistic interaction including all relevant contributions 
(and in particular hard-core repulsion and tensor) and 
2) An ab-initio computational method is employed. The 
present disagreement among different estimates of the en- 
ergy gap presently available, and which is shown in Fig. 
[U can be explained by the lack of at least one of these 
two conditions. 




FIG. 1: (color online) The 1 So pairing gap of neutron matter 
as a function of the Fermi momentum /cf computed with dif- 
ferent methods. In the figure we display works of Wambach et 
al.0], Chen et al.@|, Schulze et al.Q, Schwenk et al.QU, Cao 
et al.pH]. Gezerlis and Carlson [13] and Margueron et al.[T3|. 
All the results are compared with a BCS calculation (dashed 
line) . 



The inclusion of polarization effects shows some gen- 
eral trends in the behavior of the BCS gap. It has 
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been pointed out that the screening by the medium 
could strongly reduce the pairing strength in the 1 Sq 
channel. Actually, the diverse methods used in con- 
nection with realistic NN interactions all give a max- 
imum for A between 0.7 and 1.0 MeV at a density 
corresponding to a Fermi momentum between 0.7 and 
1.0 fm -1 ^, [nj. Recent Brueckner theory calculations 
11 1 and Hartree Fock calculations [HI give a remarkably 
larger value (~ 2.5 MeV). A recent quantum Monte Carlo 
(QMC) calculation (Til, using the AV18 NN potential pro- 
jected in the 1 Sq channel, gave in the very low density 
range corresponding to fci? <0.55 fm _1 a pairing gap sen- 
sibly larger with respect to previous results, yet lower 
than the Brueckner calculation of Ref. [ll[ , and Hartree 
Fock calculations of Ref. [l3| . 

In this letter we propose a systematic computation 
of the 1 Sq pairing gap in neutron matter as a func- 
tion of the density of the system in which both the 
conditions of using a complete realistic interaction (an 
Argonne-Urbana class potential) and an accurate ab- 
initio method (the Auxiliary Field Diffusion Monte Carlo 
within the fixed phase approximation) are fullfilled, with 
the aim of benchmarking existing results. QMC tech- 
niques have the advantage of accurately solving for the 
many-body ground state, and provide a powerful tool 
to study a wide range of systems both finite and ho- 
mogeneous. The Auxiliary Field Diffusion Monte Carlo 
(AFDMC) is particularly well suited to deal with large 
nucleonic systems 14j. In particular, the AFDMC gives 
accurate results for the properties of nucle illal , symmet- 
ric nuclear matter [ljj] and neutron matter [lj, for which 
studies performed with Green's function Monte Carlo are 
limited to 14 neutrons in a periodic box 18]. The ef- 
ficiency of AFDMC lies in the fact that spin states of 
nucleons are sampled, rather then explicitly summed, by 
means of a Hubbard-Stratonovich transformation. It can 
be extended to systems with over a hundred nucleons. 
Larger systems are needed to rule out finite size effects 
in the estimates. Moreover, we were recently able to 
overcome some technical issues related to the approxi- 
mation used to cope with the Fermion sign problem by 
turning from the constrained path approximation to a 
fixed phase approximation. This change resolved sev- 
eral issues related to discrepancies with existing Green's 
Function Monte Carlo calculations 15, TfJ lj|, 2G|. One of 
the consequences of this change is a major improvement 
in the accuracy of the AFDMC calculations for low den- 
sity neutron matter and a corresponding improvement in 
the accuracy of the calculated energy gaps over previous 
AFDMC results [21 1 . 

Our calculations for bulk neutron matter are based on 
a non-relativistic Hamiltonian of N neutrons in a periodic 
box: 
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where i and j are neutron indices, and and V%jk 
are respectively the two- and the three-nucleon interac- 
tions (namely the Argonne v' s (AV8') and the Urbana IX 
(UIX). The AV8'[22l] interactions is a simpler form of the 
more accurate Argonne vig (AV18)[23|; it contains only 
8 operators instead of 18, it preserves the same isoscalar 
part of AV18 in S and P partial waves as well as in the 
3 D± wave and its coupling to Si, and it correctly gives 
the experimental deuteron energy. The advantage of us- 
ing AV8' rather then AV18 is that it has a more suitable 
form for the AFDMC calculation. In the low density 
regime where p < 0.04 fm -3 the AV8' gives the same en- 
ergy as AV18 within 3% [20]. The Urbana UIX was fitted 
to correct the overbinding of AV18 in the ground state of 
light nuclei and to reproduce the em piri cal value of the 
equilibrium density of nuclear matter [241] . 

We consider the full nuclear Hamiltonian, instead of 
projecting in the pairing channel only, in order to include 
all the many-body correlations in the system that the 
effective bare interactions eventually miss. 

The AFDMC calculations start from a Jastrow-BCS 
trial wave function of the form 



^T = [Hfj(r ij )}$BCs(R.,S) 



(2) 



where R = {ri...rjv} and S = {si...sn} are the space 
and spin coordinates of the neutrons. The factor fj is the 
central component of the Jastrow correlation computed 
in the FHNC/SOC scheme. Its only role is to avoid the 
overlap between neutrons, and its detailed form has no 
influence on the final result. The $bcs antisymmetric 
function is built as a Pfaffian[25j] of both paired and single 
particle orbitals. Paired orbitals are defined by 
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where x is a spin function coupling two neutrons in 
the singlet state. The single particle orbitals are plane 
waves fitting Born-von Karman periodic boundary con- 
ditions. The coefficients u and v entering the paired 
orbitals are provided by a Correlated Basis Functions 
(CBF) calculation 0. In the case of even N, no single 
particle wave functions are considered in the Pfaffian, 
while if N is odd the single particle plane wave accom- 
modating the unpaired neutron is chosen in order to min- 
imize the energy of the system. Finite-size effects due to 
the truncation of the potential are reduced following the 
common procedures described in 17 1. 

Because AFDMC projects out the lowest energy state 
with the same symmetry and phase of the trial wave func- 
tion from which the projection is started, once the char- 
acter of the initial state (BCS or normal Fermi liquid) 
is given, the computed energy will refer to that partic- 
ular phase. It is therefore possible to compare the two 
equations of state, and discuss the relative stability. In 
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FIG. 2: (color online) The EOS of neutron matter in the low- 
density regime. The two calculations were performed using 
different trial wave functions modeling a normal and a BCS 
state. The fit is a guide to the eye. 



Fig. [5] we display the resulting values as a function of 
Uf normalized to the corresponding Fermi gas energy. In 
almost all of the range of densities considered the BCS 
phase is stable with respect to the normal Fermi liquid, 
although the relative energy difference never exceeds 4%. 

When the Fermi wavevector increases beyond kp=0-6 
fm -1 the normal state is energetically more favorite than 
the BCS state, with an energy difference smaller than 1% 
of the total energy. We can therefore conclude that in the 
low-density regime neutron matter is in a 1 So superfluid 
phase. In this regime the neutron-neutron interaction is 
dominated by this channel, having a scattering length of 
about a=-18.5 fm. Quantum Monte Carlo calculations 
of the EOS of dilute cold Fcrmions showed that in the 
unitary limit (when akp — > — oo), the ratio between the 
energy of the system and the energy of the Fermi gas is 
£=0.42(1) [H [2(| [13. The deviation from this asymp- 
totic value is a measure of the relevance of the details 
of the interaction in determining the equation of state 
of the system in the range of densities considered. One 
should also consider the fact that at larger densities the 
effect of pairing in scattering channels other than the 1 5o 
becomes important. 

In a full many-body calculation the superfluid gap can 
be evaluated by using the difference: 

A(N)=E(N)--[E(N + 1)+E(N-1)]. (4) 

It should be noted that the above expression is valid only 
if E(N), E(N+1) and E(N-l) are computed by keeping 
the volume V of the system fixed. This means that the 
density would be different in the N, N+l, N—l neutrons 
systems. Because our simulations are usually performed 
at fixed density, we checked the dependence of the en- 
ergy on the constraint used. Considering a number of 
particles around N = 14, which is the lowest number of 



neutrons used in the simulations, and therefore the worst 
case scenario, we evaluated the gap at fixed volume first, 
and then at fixed density. The difference in the results is 
well within statistical errors. 

Several simulations at different values of N were per- 
formed in order to evaluate the gap and the correspond- 
ing statistical error bars. A first check concerned the 
dependence of the gap estimate on the number of neu- 
tron used in evaluating the difference in Eq.(j4]). The 
values of Fermi momentum considered for the check were 
0.4, 0.6, and 0.8 fm -1 , and the numbers of neutrons were 
taken in the ranges iV=12-j-18 and N=62-^68. For each 
case, we evaluated the gap around the odd N according 
to Eq. 3J At each density A (66) (the averaged gap be- 
tween N—62 and 68) is always smaller than A(14). The 
same behavior was also observed in the QMC calculation 
of Gezerlis and Carlson [12j using the simple interaction 
projected in the pairing channel. In that paper computa- 
tions were extended also to N=92. The gap values A(66) 
and A(92) are equal within error bars and approach the 
infinite limit in the same way as in the mean field BCS 
calculation. Unfortunately in QMC simulations, in the 
absence of a correlated sampling scheme, it is impossible 
to use arbitrarily large values of N, because the gap has 
to be evaluated as the difference among total energies. 
This means that the accuracy required in the evaluation 
of the energy makes the computational time increase with 
y/~N, in addition to the N 3 standard scaling for Fermion 
simulations. 

A finite size effect might be connected to the relative 
size of the neutrons' Cooper pair and the simulation box. 
The Heisenberg uncertainty principle can be used to es- 
timate the dimension of a Cooper pair as 



where cf is the Fermi energy per particle. Taking A=2 
MeV, and fcF=0.8 fm -1 , we have therefore 5x ~ 6.6 fm, 
much smaller than the typical box size, which for 66 neu- 
trons is ~16 fm. A consequence of this analysis is also 
that relevant correlation lengths should be all contained 
in the simulation box, implying that an explicit inclusion 
in the wave function of long range effects (which are au- 
tomatically included in mean field calculations) should 
not lead to significant differences in the results. 

We report in Fig. [3] the estimate of 1 Sq superfluid 
gap as a function of the Fermi momentum fop. The 
AFDMC points are compared with results of the CBF 
calculation [5J used to determine the BCS coefficients en- 
tering the trial wave function. We also display, for sake 
of comparison, the family of more recent calculations. As 
can be seen, the AFDMC calculations give values of A 
lower than those of the CBF calculation. This behavior 
is opposite to that reported in a previous paper in which 
the same comparison was made[21|. The difference can 
be attributed both to the larger number of neutrons used 
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FIG. 3: (color online) AFDMC calculation of the So pairing 
gap of neutron matter as a function of the Fermi momentum 
hp and compared with more recent results. The AFDMC 
results are indicated by red points with the statistical error 
bars. Other results are some of those displayed in Fig. [1] and 
the blue dashed line is the CBF calculation of Fabrocini et 
al.[|. 



in the present work, and to the use of the fixed phase 
approximation instead of the constrained path approxi- 
mation to keep the sign problem under control. 

We confirm the depletion of the superfiuid gap with 
respect the BCS result. However, our results are leaning 
towards the calculations giving a maximum value of the 
gap of order 2 MeV. The other available QMC calcula- 
tions by Gezerlis and Carlson [12] differ within errorbars 
for densities corresponding to kjr <0.3fm _1 . At kp=Q.55 
fru they predict a gap about 30% smaller with respect 
to the AFDMC estimate. We believe that such difference 
comes from the fact that when increasing the density, the 
correlations induced by the interactions in channels oth- 
ers than 1 So become more and more important, and give 
a sizeable contribution to the value of the energy, and 
consequently to the gap. 
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